Effects of many-electron jumps in relaxation and conductivity of Coulomb glasses 
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A numerical study of the energy relaxation and conductivity of the Coulomb glass is presented. 
The role of many-electron transitions is studied by two complementary methods: a kinetic Monte 
Carlo algorithm and a master equation in configuration space method. A calculation of the transi- 
tion rate for two-electron transitions is presented, and the proper extension of this to multi-electron 
transitions is discussed. It is shown that two-electron transitions are important in bypassing en- 
ergy barriers which effectively block sequential one-electron transitions. The effect of two-electron 
transitions is also discussed. 

PACS numbers: 72.80.Ng 



I. INTRODUCTION 

At low temperatures, disordered systems with localized 
electrons (e. g. on dopants of compensated doped semi- 
conductors or Anderson localized states in disordered 
samples) conduct by phonon assisted hopping. The the- 
ory of this process goes back to Motti who invented the 
concept of variable range hopping. In particular he de- 
rived the Mott law for the temperature dependence of 
the conductance 

where d is the dimensionality of the sample. If Coulomb 
interactions are important one describes the system as 
a Coulomb glass due to the slow dynamics at low tem- 
peratures. As is well known (See Ref. d and references 
therein), the single particle density of states develops a 
soft gap at the Fermi level, the so called Coulomb gap. 
While this understanding of the density of states is gen- 
erally accepted, the situation is less clear when it comes 
to describing dynamics in the interacting case. Using 
the Coulomb gap density of states, which is a result of 
interactions, in the variable range hopping argument, as- 
suming that it can be used in the same way as the non- 
interacting density of states, yields the Efros-Shklovskii 
law for conductance 

While this has been observed experimentally in several 
cases it is not always the case. It remains clear that this 
is an uncontrolled approximation, and the full theoretical 
understanding of this is still missing. 

In particular, since this approach is based on a single- 
particle approximation, it neglects the possibility of cor- 
related jumps of two or more electrons. At low temper- 
atures, the system can be trapped in metastable con- 
figurations, from which it can be difficult to escape by 
single-electron transitions. By a two-electron transition 
the system can jump out of this metastable state even 
when the temperature is so low that the probability of 



making the same transition sequentially is very small be- 
cause it passes through an activated intermediate state 
with higher energy. Thus, one would expect the impor- 
tance of many-electron jumps to increase as the temper- 
ature is decreased. This was also the conclusion of some 
works^ which used a method which identifies the full 
set of low energy states of the system, and studies the 
possible transitions between them. Because the number 
of accessible states grows rapidly with increasing temper- 
ature and system size, this method is restricted to small 
systems and low temperatures. 

The importance of many-electron jumps was disputed 
by Tsigankov and Efros^, who used a kinetic Monte Carlo 
method to study the dynamics of the Coulomb glass. Us- 
ing only one-electron transitions, they confirm the Efros- 
Shklovskii law both regarding the value i in the exponent 
and also the value of Tq predicted by percolation theorj*^. 
Including two-electron transitions, they find that the two- 
electron jumps contribute about two orders of magnitude 
less to the current than the one-electron jumps. Fur- 
thermore, they find that the relative contribution of two- 
electron jumps decreases with decreasing temperature. 
They conclude that two-electron jumps are not important 
for the conductance of the Coulomb glass, contradicting 
the previous works'^'"*. We do not understand how their 
results can lead to this conclusion, since the two-electron 
jumps can be crucial in facilitating transport through 
one-electron jumps, even if their actual number is much 
less than the number of one-electron jumps. What should 
be compared is the conductance when only allowing one- 
electron jumps with the conductance when two-electron 
jumps are also included. This will be discussed in more 
detail later. 

Tsigankov and Efros^ explain their contradiction with 
the previous works"^'"* as coming from two sources: First, 
in Ref. 4 the rates of two-electron transitions were over- 
estimated because they were assumed to be independent 
of the distance between the two electrons involved in the 
transition. This assumption is not reasonable, since it 
is the Coulomb interaction which allows the remote elec- 
trons to exchange energy, and the probability of a double 
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jump should decrease with distance. Second, the method 
of identifying the full set of low energy states used in 
Refs. 0|j is numerically costly, and therefore limited to 
very small samples and low temperatures (where only 
the states with the very lowest energies are thermally ex- 
cited). Therefore the conclusions of Refs. [313 may be the 
result of small sizes and not valid for larger systems. An 
attempt to answer the second criticism was made in Ref. 
where the size dependence was studied and a scheme 
for extrapolation to infinite size was suggested. Using 
the same percolation method in configuration space as 
in Ref. (including the same expression for the many- 
electron transition rate which was questioned in Ref. 151) 
it was still found that many-electron jumps become im- 
portant at low temperatures. The difference with the 
results of Ref. |5| is explained by the fact that the method 
they used, studying all transitions between an extensive 
set of low energy states and involving up to six electrons 
moving together, was more suited to identify the cru- 
cial many-electron transitions. These could enhance the 
conductance even if they happen only very rarely since 
they could facilitate subsequent one-electron transitions. 
Since the Monte Carlo method^ becomes impracticably 
slow at low temperatures while the percolation in con- 
figuration spaced only can be applied for small systems 
and low temperatures, several questions remained unan- 
swered: Is the difference in the expression for the many- 
electron transition rate the reason for the difference in 
results? Are the two methods equivalent, or does one of 
them contain systematic errors? At what temperatures 
are the many-electron transitions important? 

In this work we try to answer these questions. Previ- 
ous works have focused on the influence of multi-electron 
transitions on conductance since this is the commonly 
measured quantity. However, at low temperatures it can 
be difficult to numerically find the conductance for two 
reasons: First, the system should be equilibrated which 
is a slow process at low temperatures. Second, to be in 
the linear response regime one must use a small potential 
difference across the sample, and the resulting current is 
so weak that one needs a long sampling time to get ac- 
curate values. Therefore we have chosen to study the 
relaxation of energy instead, this being a quantity eas- 
ily accessible in simulations. It is also well known that 
experimentally^'- the conductance is also slowly relax- 
ing, so that relaxation may be as important and relevant 
as equilibrium conductance (Although one should keep 
in mind that the experiments of Ref. H are on granu- 
lar aluminium films, and it may be that this system is 
more complex than the model discussed here accounts 
for). This allows us to go to lower temperatures using 
the Monte Carlo method, and thereby bridge the gap to 
the configuration space method. Here we report on the 
following: We give a direct calculation of the transition 
rate for two-electron jumps, to replace the unphysical one 
used in Refs. 01^ and the approximate one suggested in 
Ref. [^. A preliminary report of this part of our work 
was already presented in Ref. 9. We also discuss the ex- 



tension of this result to many-electron transitions (Sec. 
HI)) . We have numerically studied the relaxation of the to- 
tal energy of the system, comparing the evolution when 
only allowing one-electron jumps with the one where two- 
and three-electron jumps are included. We have used 
both the Monte Carlo algorithm suggested by Tsigankov 
and Efros^ and the configuration space method of Refs. 

but instead of using the percolation method we have 
studied the master equation on the set of low energy 
states, thereby eliminating any doubt on the accuracy 
of the percolation method. The details of the models 
used and the numerical procedures are given in Sec. IIIII 
The results on energy relaxation are presented in Sec. |IV] 
and some results on conductance in Sec. El 

II. MANY-ELECTRON TRANSITION RATES 

We start from the standard Coulomb gap Hamiltonian 
with a perturbation term due to tunneling 

i i<j i<j 

describing localized electrons interacting through 
Coulomb forces. c| and Ci are operators creating and 
annihilating an electron on site i, (f>i is the intrinsic 
energy of site i, which we assume to be a random variable 
uniformly distributed in the interval [—W/2,W/2], and 
Vij is the Coulomb energy. The tunneling 
amplitude tij = /oe^^''"^/" depends exponentially on the 
distance r^j, and a is the localization radius and the 
prefactor is Iq — c^/ko with k the dielectric constant. 

We consider phonon assisted tunneling due to the 
electron-phonon interaction: 

^o-ph = E E (e"''""'7q^ + h-c.) (2) 
q i 

where 6q is the phonon annihilation operator and 7q is a 
numerical factor depending on the exact phonon interac- 
tion. 

The one-electron transition rate from site i to site j is 
well known, sec for example Ref. d, 

r,, oc |7,pAr(Ai?)e-2'-../« 

where rij is the distance between the sites and AE is the 
change in energy. N{E) — l/(e^/"^ — 1) is the equilib- 
rium phonon density, for emmision processes this has to 
be replaced by N{E) + 1. We set A:^ = 1 so that temper- 
atures and energies are measured in the same units. In 
this work, following Ref. d we use the formula 

ry=ro-ie-2'->./-min (e^^^/^, l) 

where tq contains material dependent factors and energy 
dependent factors, which we approximate by their aver- 
age value; we consider it as constant and its value, of the 
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order of 10""'^^ s is chosen as our unit of time (Note that 
in Ref. a different formula was used, we do not be- 
heve that the difference is of great significance, ahhough 
it may change numerical values). 



A. Two electron transition rates 

Let us concentrate for the moment in transitions of two 
electrons on four sites and follow the method described in 
Ref. [^, based on the locator expansion in configuration 
space. We can restrict ourselves to the four sites involved 
in the transition and include the Coulomb interaction 
with the rest of the system in the site energy (pi. The 
zero-order (in the tunneling perturbation) configurations 
of two electrons on four sites are described by the states 

loo) lie) lo°> l^:> lol) (3) 

where filled circles represent sites with electrons and 
empty circles empty sites. The sites are numbered as 
III). We calculate the initial and the final states to sec- 
ond order in the tunneling, and we denote them by \ '%) 
to |°°). 

We consider phonon assisted tunneling^from the initial 
perturbed state \ %%) to the final state \ °%). The calcula- 
tion is a direct generalization of the one found in Ref. [H 
for one-electron transitions. We assume that the electron 
on an impurity is described by a hydrogen-like wavefunc- 
tion with a localization radius a and that go <C 1 where 
q is the wave vector of the phonon, so that 



vector. Further, 
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where \i) describes an electron on impurity i. After some 
algebra we find 
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Ea refers to the energy of configuration |a). This ex- 
pression corresponds to electrons jumping from site 1 to 
3 and 2 to 4. We have to add a similar expression involv- 
ing the jump from 1 to 4 and from 2 to 3. If the sites are 
at random positions, the jump with the minimum total 
hopping distance will dominate and we can neglect the 
other jumping possibilities, but if the sites are on a lat- 
tice we have to keep all the terms (and their cross terms) 
in the calculation of |Afqp. 

We further assume that qr^ ^ 1 (this may fail at 
sufficiently low temperatures) , which allows us to replace 
factors of the form cos qr^ , appearing in matrix elements 
of wavefunctions of different sites, by when integrating 
over the directions of q. Let us concentrate on the dif- 
ferent energy factors appearing in Eq. ([5]) . We first note 
that Eo c — Em m = AE is the total energy difference and 
will be equal to the energy of the phonon that is emit- 
ted or absorbed which then determines the phonon wave 
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V23 = Vl3,24 

, ; , . (6) 

is independent of the site energies. It only depends on 
the geometrical disposition of the jumps and if the sep- 
aration between sites of different jumps is much larger 
than both jumping distances it corresponds to the dipole- 
dipole interaction. The energy denominators in Eq. ([S]) 
involve the intermediate states, and it is very CPU time 
consuming to calculate these terms in numerical simula- 
tions. The divergence of these terms at certain points 
refiect the limitation of the perturbation theory rather 
than any physical effect. Therefore we want to cut off 
this divergence and replace the fraction by 1 when it is 
larger than 1. The energy differences are of the order 
of the disorder W. Therefore the terms in the brackets 
are also never very small. Since they are also tempera- 
ture independent we propose to set these terms equal to 
1/W. We believe that this is of no physical consequence, 
and will not affect the results qualitatively. Taking into 
account these approximations and using Fermi's golden 
rule we arrive at the expression for the transition rate 
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(7) 

where tq is the same unit of time as for one-electron tran- 
sitions and we assume that the average energy difference 
of intermediate states is W/2. More details on the deriva- 
tion of this equation were given in Ref. 



B. Three and more electron transition rates 

In this case one has to assume that the interaction is 
weak and do perturbation theory in both the interaction 
and the hopping term. As for two-electron jumps, it is an 
excellent approximation for sites at random to consider 
only the transitions with the minimum total hopping dis- 
tance Ri^j. There are many of them, differing in the or- 
der of the one electron moves and on the jump directly 
excited by the phonon. The final expression for the tran- 
sition rate is very complex, but its most important factor 
is easy to get^^. It is the probability of finding a phonon 
of energy AE times the product of the overlap integrals 
of the hops of all electrons involved in the transition. The 
many electron transition rate can then be approximated 
by 



,-AE/T 



(8) 



7 is a measure of the importance of the interaction energy 
compared to the disorder energy and n is the number of 
electrons participating in the process. 

The use of equation ([8]) for the transition rates in nu- 
merical simulations may overestimate the importance of 
correlated hops since it doublecounts the effects of exci- 
tations well separated one from each other. Although 
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one-electron excitations should dominate in this case, 
since many-electron excitations should only be important 
when single excitations have positive energies, while the 
combined excitation has negative energy, it is convenient 
to get rid of this problem. We can do that by substi- 
tuting the constant 7 by a prefactor similar to the one 
obtained for two-electron transitions, equation ([7]). It is 
difficult to get a closed expression for this prefactor, and 
we propose an empirical approach that it is practical for 
numerical purposes and that we think incorporates the 
relevant physics of the problem. A requirement for this 
prefactor is that it should vanish when one of the transi- 
tions is very far from the others. A suggestion satisfying 
this requirement is the sum of all the products of rt — 1 
different interaction energies between any pair of single 
electron transitions, like Vl3.24 in (H]). Each of these terms 
must be divided by a factor proportional to the disorder 
energy as in the two-electron case. This proposal corre- 
sponds to exciting one of the hops by the phonon and 
the rest by the dipole-dipole interaction in all the possi- 
ble ways. For three electron transitions we take as the 
preexponential 

-^022 22 22 

■1^77(^14, 25^25, 36 + ^25,36^14,36 + ^14,25^14,36)- i^) 



III. MODEL FOR NUMERICAL SIMULATIONS 

We use the standard tight-binding Coulomb gap 
Hamiltonian^: 

(ni — K)(nj — K) , ^, 

H = Y.e,n, + Y.— ^ (10) 

■ ^ - ' ii 

K being the compensation. We take the number of elec- 
trons to be half the number of sites. The sites are ar- 
ranged in two dimensions both on a lattice and at ran- 
dom, but in the latter case with a minimum separa- 
tion between them, which we choose to be O.OSr where 
= L"^ /N . We implement cyclic boundary conditions 
in both directions. We take e^/r as our unit of energy 
and r as our unit of distance. 



A. Monte Carlo algorithm for lattice systems 

For the two-electron transition rate we use Eq. ([7]) for 
sites at random and the extension that includes the two 
jumping possibilities when sites are on a lattice. As for 
one-electron transitions, the rate is split in one energy 
dependent (or activation) term, F"^ and one distance de- 
pendent (or tunneling) term F"^. This means that we can 
use the hybrid algorithm of Tsigankov and Efros'^. 

The program first calculates and stores the distance 
dependent part of the rates. For the one electron jumps, 
the tunneling parts of the rate for all jumps in the square 
|Aa;|, jAyl < Li some maximal jump length (in the nu- 
merics Li = 10 lattice units) are calculated. The sum 



pTotai _ Y^Yt,! is also stored. For the two electron 
jumps, all coordinates are relative to the initial position 
of the first electron. The following algorithm is used: 

• The final position of the first electron is selected 
in the square |Aa;|, |A?/| < L2 where the size L2 
can be reasonably chosen to be about half of ii 
since the distances each electron jumps are added 
together to find the rate (In the numerics L2 = i 
lattice units). 

• The initial position of the second electron is se- 
lected in the square |a;2|, I2/2I < D2 (In the numer- 
ics D2 — 5 lattice units) . The initial position of the 
second electron can not be either the initial or the 
final position of the first electron. 

• The final position of the second electron is selected 
in the square I Aa;2 1, |Ay2| < -^2- The final position 
of the second electron can not be the initial or final 
position of the first electron. 

• The tunneling part of the rate for this transition is 
calculated according to the formula 

Ft,2 = -£'1^13,24 + ^2^14,23 + -£'1^2143,24^44,23 

where 

J^^ = g-2(r-i3-|-r24)/o^ ^2 _ g-2(ri4-|-T-23)/a 

• The rates for all these transitions are stored, and 
the sum of aU F|^°2*''' = J2 ^t,2 is calculated. 

When the Monte Carlo algorithm is running, it will do 
the following steps in order to select which transition to 
make. 

• It is decided whether to attempt a one electron 
transition (probability F^°*''7(r?°i*''' + r|;°*°0) 
or a two electron transition (probability 

-nTotal //-nTotal , -nTotal\\ 
T,2 / \^ T,l "T T,2 )) 

• If it is a one electron transition we follow the usual 
procedure of Tsigankov and Efros^. One occupied 
site is selected randomly. Then the final site is se- 
lected at random but with weights given by the 
probability Ft,i/F|;°*'''. If the final site IS occu- 
pied, the transition is rejected. If it is empty it is 
accepted with probability min(l, e^'^^/^). 

• If it is a two electron transition an occupied site 
is selected at random. Then a certain two elec- 
tron jump is selected weighted by the probability 

Tt,2/^tT'^- If ^^^^ °f electron 
is occupied, the initial site of the second electron 
is empty or the final site of the second electron is 
occupied, the transition is rejected. If not, the tran- 
sition is accepted with probability min (l,e-'^^/^). 
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B. Master equation method for sites at random IV. RESULTS ON ENERGY RELAXATION 



We used a numerical algorithm to obtain the ground 
state and ~ 10^ lowest energy many-particle configura- 
tions of the system. The algorithm is an improved ver- 
sion of the algorithm in Ref. 6 and it was described in 
detail in [Til It consists of the following two stages. In 
the first stage, we repeatedly start from states chosen at 
random and relax each sample by means of a local search 
procedure. In an iterative process, we look for neighbors 
of lower energies and always accept the first such state 
found. The procedure stops when no lower energy neigh- 
boring states exist, which insures stability with respect to 
all one-electron jumps and compact two-electron jumps. 
We then consider a set of metastable sates found by 
the process just described and look for the sites which 
present the same occupation in all of them. These sites 
are assumed to be frozen, i.e. they are not allowed to 
change occupation, and the relaxation algorithm is now 
applied to the unfrozen sites. The whole procedure is 
repeated until no new frozen sites are found with the set 
of metastable states considered. 

In the second stage, we complete the set of low-energy 
configurations by generating all the states that differ by 
one- or two-electron transitions from any configuration 
stored. In order to speed up this process, which is very 
CPU time consuming, we again assume frozen and un- 
frozen sites and in first place look for neighboring con- 
figurations by changing the occupation of unfrozen sites 
only. We later relax this restriction in the final stage. 

We consider 20 different realizations of a system with 
500 sites. We obtain the 200 000 lowest energy configu- 
rations for each realization. We then obtain all one- two- 
and three electron transitions between all these configu- 
rations, establishing a dynamical matrix that evolves the 
system in time according to this master equation. We 
choose for the initial state the mixture of all configura- 
tions with weights equal to the Boltzmann factor for a 
temperature 40 times larger than the real one. 

We have developed a renormalization procedure to 
propagate in time efficiently. It takes advantage of the 
fact that the distribution of relaxation times is exponen- 
tially broad and at large times the short time processes 
must have already equilibrated some sets of configura- 
tions. At a given time, we calculate the current through 
each transition and if it is very small, relative to the tran- 
sition rate and the occupation probability, we assume 
that the two configurations involved are in thermal equi- 
librium and can consider then as part of a cluster. We 
recalculate the transition rates between this new clus- 
ter and the rest of the system. The time step used in 
the numerical propagation is calculated dynamically and 
increases drastically as we form more clusters or larger 
clusters. 



A. Relaxation using Monte Carlo on the lattice 
model 

To see the effect of two electron transitions on re- 
laxation we did the following. For one sample of size 
100 X 100 and = 2 we ran relaxation from an initial 
random state at three different temperatures, T = 0.01, 
0.005 and 0.001. For each temperature we ran from the 
same initial state 10 (20 for T = 0.001) different time 
evolutions (different random seeds for the Monte Carlo 
evolution). For each case we ran the simulation both 
allowing only one electron transitions and including one 
and two electron transitions. The results are shown in 
Fig. [TJ It is clear that as the temperature decreases, the 
difference between the relaxation rates in the one- and 
two-electron cases increases. To confirm that the results 
are general and the sample sufficiently large we did one 
set of 10 time evolutions on the same sample but start- 
ing from a different initial state and one set on a different 
sample. The same behaviour was seen in all cases. 

To see more clearly the importance of the two electron 
jumps we can look at one particular relaxation graph and 
mark the points where two electron jumps occur (Figure 
O the two graphs, Figure [2^) and Figure [5}d) correspond 
to two different Monte Carlo evolution of the same sam- 
ple and initial state). In the figure the curve is the en- 
ergy, while the points mark the time when a two electron 
jump was performed. The red points represent the final 
energy after the transition, while the green points rep- 
resent the energy of the intermediate state if this jump 
was to be replaced by sequential one electron jumps. The 
temperature was T = 0.001, and as we can see the in- 
crease in energy to the intermediate state is sometimes 
more than two orders of magnitude larger that this. This 
means that the probability of this process occurring se- 
quentially is extremely small. As can be seen from the 
figure, there are clear correlations between the occurrence 
of two electron jumps and steps in the relaxation graph. 
This means that the two electron jumps are essential in 
overcoming barriers in the relaxation path and give a 
contribution to the relaxation rate even if the number of 
two electron jumps can be a small fraction of the total 
number of jumps. Sometimes (at long times in Fig ^jp) 
there can occur soft two electron excitations which are 
then jumping back and forth between the two configura- 
tions like soft dipoles in the one electron case. These give 
large contributions if we try to count the number of two 
electron transitions, but are not important for relaxation. 

To better measure the importance of the two electron 
jumps on the relaxation, we do the following. For each 
decade in time we see how much the energy was reduced 
in the two electron jumps (or by one electron jumps im- 
mediately following a two electron jump) and compare 
this to the total relaxation of energy during this time. 
We then get Fig. [3] and as we can see, the two electron 
jumps increase in their importance for relaxation, and 
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FIG. 1: (Color online) Energy relaxation as function of time at a) T = 0.001, b) T = 0.005 and c) T = 0.01. Averages are 
shown for one and two electron transitions, error bars are standard deviation of the mean. 
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FIG. 2: (Color online) The correlation between the points 
where two electron jumps occur and steps in the energy re- 
laxation graph. The blue curve is the total energy of the 
system. The red points (-I-) are the final energies after two 
electron jumps and the green points (*) the energy of the in- 
termediate state if this jump was to be replaced by sequential 
one electron jumps. T — 0.001. The two graphs, a) and b), 
correspond to two different Monte Carlo evolution of the same 
sample and initial state 
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FIG. 3: The fraction of relaxation due to two electron tran- 
sitions. The figure is an average over 6 different Monte Carlo 
evolutions on one sample at T = 0.001. 



at the later times in the simulation, they are responsible 
for most of the relaxation. At long times the fraction 
slightly exceeds 1 which is due to the fact that the en- 
ergy of the system was not stored after each jump so the 
energy reduction in a two electron jump can be slightly 
overestimated if it was preceeded by one electron jumps 
which increased the energy. 



B. Relaxation using master equation on the 
random sites model 



We have calculated the average energy, with respect 
to the ground state energy, as a function of time. In 
Figure |3I we plot the results for one-electron and two- 
electron relaxations at a temperature T — 0.002 for a 
system with 500 sites. We consider this small size in 
order to have configurations extending over a relatively 
large energy range. The dotted line is the result when 
only one-electron relaxation is considered. The continu- 
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FIG. 4: (Color online) Energy relaxation as function of time 
for one-electron jumps (black dotted curve), one- and two- 
electron jumps without any prefactor (green dashed curve) 
and with the prefactor of Eq. (blue continuous curve). 



FIG. 5: (Color online) Energy relaxation as function of time 
for one-electron jumps (black dotted curve), one- and two- 
electron jumps (blue dashed curve) and up to three-electron 
jumps (red continuous curve). 



ous and the dashed curves correspond to relaxation by 
one- and two-electron jumps. In the dashed case we have 
not included any prefactor in the expression for the tran- 
sition rates, while in the continuous curve we use Eq. 
([7|). We first note how relaxation by one-electron jumps 
alone is far from complete. The system gets easily stuck 
in metastable states even for the relatively small system 
size considered. The inclusion of two-electron jumps is 
almost negligible at short times, where one can always 
lower the energy by one-electron jumps. But at larger 
times, two-electron transitions are really needed to over- 
come the energy barriers. 

We also note that if we do not include the right pref- 
actor in the two-electron rate we are overestimating their 
effects, specially at short relaxation times, because we are 
double counting some excitations. We have checked that 
the results for two-electron contributions do not change 
if we only include those transitions with negative inter- 
action energy. This result in an important reduction in 
the number of many-electron excitations to include in the 
simulations. In future calculations of many-electron ef- 
fects it will be convenient to take advantage of this result 
and to explore more drastic reductions in the number of 
relevant excitations. 

At a time of roughly IO^tq we have already practi- 
cally reached the thermal equilibrium when up to two- 
electron transitions are considered. We expect this time 
to increase drastically with system size. In figure [5] we 
have represented energy relaxation by one-electron jumps 
(dotted curve) , by up to two-electron hops (dashed curve) 
and by up to three-electron jumps (continuous curve). 
For three-electron jumps we have used for the prefactor 
the sum of all the different products of dipole-dipole con- 
tributions. We note that the inclusion of three-electron 
jumps does not affect much the results for the size consid- 



ered. We expect that their effects will be more important 
for larger sizes, which will contain larger energy barriers 
and a more complex energy landscape. As we include 
transitions of more particles, the energy relaxation curve 
seems to approach a logarithmic behavior. 



V. RESULTS ON CONDUCTIVITY 

We also studied conductivity, comparing the cases with 
and without two electron jumps. At temperatures T > 
0.04 we ran four different samples of size 100x100. For 
T < 0.04 we used four samples of size 200x200 (this be- 
cause we know that at low temperatures we see finite size 
effects in the conductance up to L — 100). The electric 
field was T/10, which should be in the ohmic regime. In 
each case we ran the simulation for 10^ accepted jumps 
and checked that this was sufficient to obtain a straight 
line of transferred charge as function of time. The con- 
ductance is then given by the slope of this line. The re- 
sults are shown in Fig. [51 We see that we reproduce the 
Efros-Shklovskii law for the conductivity and that there 
is no significant difference when including two-electron 
transitions. 



VI. DISCUSSION 

Comparing our results with the previous works^^'' it 
seems that all fall into the same coherent picture. By 
focusing on relaxation, we were able to apply the Monte 
Carlo method to temperatures comparable to the ones 
were the master equation approach can be used. Then 
we see similar behavior in both the cases. The energy 
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FIG. 6: (Color online) Conductance as function of temper- 
ature showing the Efros Shklovskii law. Blue points (*) are 
with one electron jumps, red points (+) include also two elec- 
tron jumps. Four samples are shown for each temperature, 
the spread indicates the uncertainty. The lines are the aver- 
ages. 



relaxes faster when two electron transitions are included. 
A more direct comparison of the two methods is difficult 
since the models differ in several respects. In the Monte 
Carlo method we prefer to use a lattice since it reduces 
the computational effort while in the master equation we 
are resticted to small samples and prefer a random site 
model since we believe that lattice effects are more severe 
for small systems. Also, the initial states are different 
in the two cases, since in the master equation approach 
we need to take as initial state some combination of the 
states in the low-energy set we are working on. These are 
already very low-energy states, and different in structure 
from the random states used in the Monte Carlo method. 
However, we find that our results convincingly show that 
both methods give similar results at low temperatures. 



and that there are no systematic errors which affect one 
or the other method. Furthermore, if we look at Figure |3] 
and compare the two graphs including two electron tran- 
sitions, with and without the prefactor in Eq. ([T]), we find 
that although the omission of the prefactor overestimates 
the importance of two electron jumps the results remain 
qualitatively the same. Thus, the concern of Tsigankov 
and Efros^ that this overestimation changes the results 
qualitatively seems unfounded. We may then still be- 
lieve the results of Ref d, at least on the qualitative level. 
Comparing our Figure [S] with Figure 2 of Ref. d, both 
calculations of the conductivity using the same Monte 
Carlo method, we find that the two agree closely (a de- 
tailed comparison shows a small shift in the values but 
the slope of the line remains the same) . Figure 4 of Ref. 
|6| gives the corresponding results using the configuration 
space approach. We see that although the configuration 
space method finds a difference in the conductivity when 
including two electron jumps, this difference is small, at 
the level of the statistical error, for I/a/T < 8 which is 
where we have results using the Monte Carlo method. 
If we compare with our Figure (TJ we see that at these 
temperatures we can not see any significant effect of two 
electron jumps on relaxation either. We therefore con- 
clude that two electron jumps will only be important at 
lower temperatures, and we believe that we would also 
see this in Monte Carlo simulations if these could be per- 
formed at sufficiently low temperature. 
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